#Required libraries

library(tidyverse)
library(broom)
library(lubridate)
library(viridis)
library(RColorBrewer)

#Data import

abn <- read.csv("./Data/Abundance_280820.csv")
transects <- read.csv("./Data/Transects.csv")

# adding relevant variables

abn <- abn%>%
  mutate(Year = year(dmy(Date)))%>% #adding year of sampling
  mutate(Study = case_when(Year %in% c(1998, 1999, 2000, 2001) ~ "1998 - 2000", #previous study
                           Year %in% c(2016, 2017) ~ "2016 - 2017"))%>% #current study
  mutate(Transect = as.character(Transect),
         Transect = ifelse(Transect == "T6A" | Transect == "T6B", "T6", Transect),
         Transect = as.factor(Transect))

transects <- transects%>%
  mutate(Transect = as.factor(Transect))

# cleaning abundance
abn <- abn%>%
  filter(!is.na(Abundance) &
         Specific.Epithet != "" & #removing samples with no species
         Specific.Epithet != " " &
          Specific.Epithet != "Unidentified" & #removing unidentified species
         substr(Specific.Epithet, 1, 1) != "?")%>% #83 obs removed
  droplevels()

theme_set(theme_bw()+
            theme(legend.position = "top",
                  text = element_text(size = 20)
                  )
          )#setting ggplot theme

# function for t-test reporting

tstats <- function(x){
  x%>%
    rename('1998-2001' = estimate1,
         '2016-2017' = estimate2)%>%
    dplyr::select(`1998-2001`, `2016-2017`, parameter, statistic, p.value)
}

abn%>%
  group_by(Study, Taxa, Transect)%>%
  summarise(n = length(unique(Date)))%>%
  spread(Transect, n)

abn%>%
  group_by(Study, Taxa, Transect)%>%
  summarise(n = sum(Abundance))%>%
  spread(Transect, n)

# sample info

samp_birds <- abn%>%
  left_join(transects, by = c("Transect"))%>%
  filter(Taxa == "Birds")%>%
  dplyr::select(Study, Taxa, Transect, Habitat)%>%
  distinct()%>%
  mutate(sample = paste(Study, "_", Transect))

samp_butterfly <- abn%>%
  left_join(transects, by = c("Transect"))%>%
  filter(Taxa == "Butterflies")%>%
  dplyr::select(Study, Taxa, Transect, Habitat)%>%
  distinct()%>%
  mutate(sample = paste(Study, "_", Transect))

# Similarity matrix by studies

sim_birds <- abn%>%
  filter(Taxa == "Birds")%>%
  group_by(Study, Specific.Epithet)%>%
  summarise(n = sum(Abundance, na.rm = T))%>%
  mutate(sample = Study)%>%
  ungroup()%>%
  dplyr::select(sample, Specific.Epithet, n)%>%
  spread(Specific.Epithet, n, fill = c(n = 0))%>%
  column_to_rownames(var = "sample")

sim_butterfly <- abn%>%
  filter(Taxa == "Butterflies")%>%
  group_by(Study, Specific.Epithet)%>%
  summarise(n = sum(Abundance, na.rm = T))%>%
  mutate(sample = Study)%>%
  ungroup()%>%
  dplyr::select(sample, Specific.Epithet, n)%>%
  spread(Specific.Epithet, n, fill = c(n = 0))%>%
  column_to_rownames(var = "sample")

# bray - curtis

data_frame(
  Taxa = c("Birds", "Butterflies"),
  Dissimilarity = c(as.numeric(vegdist(sim_birds, method = "bray")),
                    as.numeric(vegdist(sim_butterfly, method = "bray")))
)

# ordination

NMDS_birds <- metaMDS(veg_birds, distance = "bray", k = 2, model = "local")

NMDS_butterfly <- metaMDS(veg_butterfly, distance = "bray", trymax = 1000, k = 2, model = "local")


# Extracting scores and adding habitat vars

ord_birds <- data.frame(scores(NMDS_birds))%>%
  rownames_to_column("sample")%>%
  left_join(samp_birds, "sample")

ord_butterfly <- data.frame(scores(NMDS_butterfly))%>%
  rownames_to_column("sample")%>%
  left_join(samp_butterfly, "sample")

ord <- bind_rows(ord_birds, ord_butterfly)

adonis(veg_birds ~ Study, data = samp_birds)

adonis(veg_butterfly ~ Study, data = samp_butterfly)

# plotting

ord%>%
  ggplot(aes(NMDS1, NMDS2, col = Study, fill = Study))+
  stat_ellipse(level = 0.95, geom = "polygon", linetype = "dashed", size = 1, alpha = 0.3)+
  geom_point(aes(shape = Habitat), size = 3)+
  guides(shape = guide_legend(nrow = 2, byrow = T))+
  scale_color_brewer(palette = "Set1")+
  scale_fill_brewer(palette = "Set1")+
  facet_wrap(~Taxa)+
  theme(legend.box = "vertical")

ggsave(last_plot(), filename = "./Figures and Tables/figure3.png", height = 6, width = 8, dpi = 300)

season <- abn%>%
  mutate(Month = month(dmy(Date), label = T))%>%
  group_by(Study, Month, Taxa)%>%
  summarise(n = sum(Abundance, na.rm=T))%>%
  group_by(Study, Taxa)%>%
  mutate(N = sum(n),
         rel.prop = n/N)

season%>%
  ggplot(aes(Month, rel.prop, col = Taxa, shape = Taxa))+
  geom_point(size = 3)+
  geom_path(aes(group = Taxa), size = 1, linetype = "dashed")+
  scale_color_brewer(palette = "Set1")+
  facet_wrap(~Study, ncol = 1)

library(raster)
library(sf)

# Raw point data

points <- read.csv("./Data/sites.csv")%>%
  drop_na()%>%
  dplyr::select(-Location)

# Converting to SF points

points_sf <- st_as_sf(points, coords = c('Lon', 'Lat'))

st_crs(points_sf) <- 4326

# Making transect lines

lines <- points_sf%>%
  mutate(Transect = as.factor(Transect))%>%
  group_by(Transect)%>%
  summarise()%>%
  st_cast("LINESTRING")

# getting study extent

study_ext <- st_bbox(points_sf)

ext <- st_as_sfc(study_ext)

# getting osm map

library(osmdata)

# Tamhini WLS

tamhini_raw <- opq(bbox = study_ext, timeout = 100)%>%
  add_osm_feature(key = 'boundary', value = 'protected_area')%>%
  osmdata_sf()%>%
  unique_osmdata()
  
tamhini <- tamhini_raw$osm_multipolygons

st_crs(tamhini) <- "+init=epsg:4326"

# Roads

tamhini_area <- st_bbox(tamhini)

roads_raw <- opq(bbox = tamhini_area, timeout = 100)%>%
  add_osm_feature(key = 'highway')%>%
  osmdata_sf()%>%
  unique_osmdata()
  
roads <- roads_raw$osm_lines

st_crs(roads) <- 4326

# Landuse

landuse_raw <- opq(bbox = tamhini_area, timeout = 100)%>%
  add_osm_feature(key = 'natural')%>%
  osmdata_sf()%>%
  unique_osmdata()
  
forest <- landuse_raw$osm_polygons%>%
  filter(natural %in% c("wood", "tree"))

wood <- landuse_raw$osm_multipolygons%>%
  filter(natural == "wood")

water <- landuse_raw$osm_multipolygons%>%
  filter(natural == "water")

water_small <- landuse_raw$osm_polygons%>%
  filter(natural == "water")

st_crs(forest) <- 4326

# waterways

waterway_raw <- opq(bbox = tamhini_area, timeout = 100)%>%
  add_osm_feature(key = 'waterway')%>%
  osmdata_sf()%>%
  unique_osmdata()
  
waterways <- waterway_raw$osm_lines

st_crs(waterways) <- 4326

# elevation 

library(elevatr)
library(rgdal)

tam_elev <- get_elev_raster(locations = tamhini, z = 8) 

tam_contours <- rasterToContour(tam_elev, levels = seq(0, 1000, 100))

# plotting site map

library(tmap)
library(cowplot) # for plot_grid() function - good to arrange multiple plots into a grid
library(grid)

# Study site

tm_shape(forest)+ #forested lands
  tm_polygons(col = "#81AF81", alpha = 0.5, border.col = NULL)+
tm_shape(wood)+ #forested lands
  tm_polygons(col = "#81AF81", alpha = 0.5, border.col = NULL)+
tm_shape(tam_contours)+
  tm_lines(col = "white", lty = "dotted")+
  tm_text(text = "level", auto.placement = T)+
tm_shape(water)+
  tm_polygons(col = "#77A1CB", border.col = NULL)+
tm_shape(waterways)+
  tm_lines(col = "#77A1CB", lwd = 2)+
tm_shape(roads)+
  tm_lines(alpha = 0.5, lwd = 1.5)+
  tm_text("ref", size = 0.3, remove.overlap = T)+
tm_shape(tamhini)+
  tm_borders(lwd = 2, col = "#165916", lty = "dashed")+
  tm_text("name", size = 0.5)+
tm_shape(lines, is.master = T, bbox = study_ext+c(-0.005, -0.005,0.005,0.005))+
  tm_lines(col = "red", lwd = 2)+
  tm_text("Transect", size = 0.5, auto.placement = T, remove.overlap = T)+
  tm_graticules(n.x = 5, lines = F, labels.inside.frame = T)+
  tm_scale_bar(position = c("left", "top"), text.size = 0.5, size = 0.5) + 
  tm_compass(position = c('right', 'top'), size = 2)+
  tm_layout(scale = 1.5, legend.position = c("left", "bottom"), bg.color = "#ffdb9c", asp = 1.5)

fig1a <- grid.grab()

# Tamhini WLS

tm_shape(forest)+ #forested lands
  tm_polygons(col = "#81AF81", alpha = 0.5, border.col = NULL)+
tm_shape(tam_contours)+
  tm_lines(col = "white", lty = "dotted")+
  tm_text(text = "level", auto.placement = T)+
tm_shape(water)+
  tm_polygons(col = "#77A1CB", border.col = NULL)+
  tm_text(text = "name", remove.overlap = F, size = 0.5)+  
tm_shape(water_small)+
  tm_polygons(col = "#77A1CB", border.col = NULL)+
  tm_text("name", remove.overlap = F)+  
tm_shape(waterways)+
  tm_lines(col = "#77A1CB", lwd = 2)+
tm_shape(roads)+
  tm_lines(alpha = 0.5, lwd = 1.5)+
tm_shape(tamhini, is.master = T, bbox = tamhini_area+c(-0.01, -0.01,0.01,0.01))+
  tm_borders(col = "#165916", lty = "dashed")+
  tm_text("name", size = 0.5)+
tm_shape(ext)+
  tm_borders(col = "red")+
  tm_scale_bar(position = c("left", "top"), text.size = 0.5) + 
  tm_layout(scale = 1, legend.position = c("left", "bottom"), bg.color = "#ffdb9c", asp = 1)

fig1b <- grid.grab()

# legend

tm_shape(forest)+
  tm_polygons(col = "#81AF81", alpha = 0.5, border.col = NULL)+
tm_shape(ext)+
  tm_borders(col = "red")+
  tm_add_legend(type = "fill",col = c("#81AF81", "#77A1CB", "#ffdb9c"), 
                labels = c("Forest", "Water", "Barren Land"), title = "Legend")+
  tm_add_legend(type = "line", col = c("black", "#165916", "#77A1CB", "red"), 
                labels = c("Roads", "WLS boundary", "Waterway", "Transect"), 
                lwd = c(3,3,3,3),
                lty = c(1, 2, 1, 1))+
  tm_layout(legend.only = T, asp = 1, legend.position = c(0.3,0.2), scale = 2)

legend <- grid::grid.grab()

# making inset map

tamext <- st_as_sfc(tamhini_area + c(-.5,-.5, .5, .5)) 

india <- read_sf("./Data/India-States.shp")

tmap_options(check.and.fix = TRUE)

tm_shape(india)+
  tm_polygons(col = "#ffdb9c")+
tm_shape(tamext)+
  tm_borders(lwd = 1, col =  "red")+
  tm_layout(scale = 1.5, bg.color = "#77A1CB", asp = 1)

fig1c <- grid.grab()

# making final map

library(gridExtra)

fig1 <- grid.arrange(fig1a, fig1b, fig1c, legend, layout_matrix = rbind(c(1, 1, 1),
                                                          c(1, 1, 1),
                                                          c(2, 3, 4)))

ggsave(fig1, filename = "./Figures and Tables/figure1.png", height = 8, width = 8)